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Abstract 


We have analyzed the transferability of a previously proposed Buckingham 
repulsion-dispersion intermolecular potential for the explosive hexahydro-1,3,5-trinitro-1,3,5-s- 
triazine (RDX) (Sorescu, D. C., B. M. Rice, and D. L. Thompson, Journal of Physical 
Chemistry B, vol. 101, p. 798,1997) to predict the crystal structures (within the approximation 
of rigid molecules) of a database of 30 nitramines. These include acyclic, monocyclic, and 
polycyclic molecules. It is shown that the proposed potential model is able to accurately 
reproduce the crystallographic structures and lattice energies (where available) of these crystals. 
For the majority of these crystals, the best agreement with experimental structural and energetic 
data is obtained in those cases when the electrostatic charges have been determined using ab 
initio methods that include electron correlations effects (namely MP2 and B3LYP). The use of 
the electrostatic charges calculated at the Hartree-Fock level results in large deviations of the 
predicted lattice energies from the experimental values. The deviations of the lattice energies can 
be significantly decreased by scaling the electrostatic charges with a constant factor. 
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1. Introduction 


Atomistic simulation is increasingly gaining acceptance as a practical research tool in the 
investigation of the behavior of condensed-phase materials. In addition to providing information 
that is difficult or impossible to measure, prediction leads to a reduction in unnecessary 
measurement or synthesis of candidates in the course of design of new materials. However, the 
power of atomistic simulation can only be realized if the description of the molecular system is 
accurate. The development of accurate intermolecular potentials is not a simple, straightforward 
procedure. Substantial work has been directed toward determining both simple functions that 
make large-scale simulation realizable, and correct parameterization such that the physical 
properties of the materials are properly described. In this work, an intermolecular potential that 
accurately describes nitramine crystals is presented. Also, the potential parameters effect on 
predictive ability is investigated. 

In initial studies of nonreactive processes in the nitramine explosive hexahydro-1,3,5,-trintro- 
1,3,5-s-triazine (RDX), an intermolecular potential energy function was developed that would 
accurately reproduce the structure of the a-form of the RDX crystal [1]. This potential is 
composed of pairwise atom—atom (6-exp) Buckingham terms with explicit inclusion of the 
electrostatic interactions between the charges associated with the atoms of different molecules. 
The parametrization of the potential function was done such that molecular packing calculations 
(MP) reproduced the experimental structure of the crystal and its lattice energy. 
Isothermal-isobaric molecular dynamics simulations (NPT-MD) using this potential energy 
function predicted crystal structures in excellent agreement with the experimental data [1]. 

It has been shown that this interaction potential energy function is transferable to two other 
nitramine crystals: the polycyclic nitramine 2,4,6,8,10,12-hexanitrohexaazaisowurtzitane (HNIW) 
[2] and the monocyclic nitramine l,3,5,7-tetranitro-l,3,5,7-tetraazacyclooctane (HMX) [3]. Both 
MP and NPT-MD simulations predict geometrical parameters in good agreement with the 
experimental values for the different polymorphs of the HNIW and HMX crystals [2, 3]. 
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Furthermore, the calculations indicate a stability ranking for HNIW in agreement with 
experimental measurements [4]. 

The success of these potential energy parameters in describing the RDX crystal and different 
phases of the HMX and HNIW crystals provides impetus for further investigations to determine 
the limits of the transferability of this interaction potential. Toward this end, MP calculations of 
30 nitramine crystals are reported here. This set of crystals is composed of monocyclic, 
polycyclic, and acyclic nitramine molecules. There was a particular interest to see if the 
geometrical and energetic parameters for these crystals could be reproduced by this proposed 
model. 

One of the main factors that contributes to the quantitative description of the molecular 
packing in a crystal is related to the representation of the electrostatic interactions. It was shown 
more than a decade ago [5, 6] that increased accuracy in structural predictions of the molecular 
crystals and in transferability of the potential parameters can be achieved by explicit use of the 
electrostatic interactions between the charges associated with the atoms. For example, many of 
the available force fields such as Amber [7], ECEPP [8], or Dreiding [9] use these kinds of 
potential terms in models of organic, biological, and main group inorganic crystals. Further 
improvement of the description of the electrostatic forces between molecules, particularly in 
crystals with substantial anisotropies, can be achieved by using sets of point multipoles (charge, 
dipole, quadrupole, etc.) on every atomic site. This distributed-multipole representation has been 
shown to be successful in the modeling of the crystal structures of polar and hydrogen-bonded 
molecules [10]. 

The present study has found that, as in the cases of the RDX [1], HNIW [2], and HMX [3] 
crystals, the set of 30 crystals considered here can be accurately represented using the 
Buckingham potential plus Coulombic interactions. The assignment of the electrostatic charges 
poses a problem in that the atom-centered monopole charge is not an observable quantity and 
cannot be obtained directly from either experiment or first principles calculations. Currently, 
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there are several schemes for evaluation of charges by empirical partition or by using a quantum 
mechanically derived wave function [11-13]. The Coulombic terms have been determined 
through fitting of partial charges centered on each atom of the molecules to a quantum 
mechanically derived electrostatic potential [13]. An investigation has been done on how the 
geometrical and energetic parameters predicted in MP calculations depend on charges determined 
from ab initio methods, which do or do not include electron correlation effects. Specifically, 
different sets of charges derived from the Hartree-Fock (HF) wave function [14] and from 
methods that employ electron correlations such as second-order Moller-Plesset (MP2) 
[15-18] and B3LYP [19, 20] have been used. 

The studies described here represent the first stage in the development of a general model for 
nitramine crystals. The main limitations of the present model are due to the assumption of rigid 
molecules but further refinement of this model can be made to include the effects of 
intramolecular motions, particularly of low-frequency torsional motions of the nitro groups and 
the ring. 

The organization of the paper is as follows. In section 2, the intermolecular potential used to 
simulate the nitramine crystals is presented. In sections 3 and 4, the molecular packing methods 
and results, respectively, are described. The main conclusions are summarized in section 5. 

2. Intermolecular Potential 


The central problem in classical simulations of molecular crystals is the construction of 
re alis tic potentials that accurately predict the structural and thermochemical parameters. In this 
paper, the same general model is employed for the atom-atom potentials that proved to be 
successful in modeling of the RDX, HNIW, and HMX crystals [1-3]. In particular, it is assumed 
that (1) the intermolecular interactions depend only on the interatomic distances; (2) the 
interaction potential can be separated in contributions identified as van der Waals and 
electrostatic, and (3) the same type of van der Waals potential is used for the same type of atoms, 
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independent of their valence state. Moreover, in the present case, the transferability of the 
potential parameters determined for the case of RDX crystal to all the nitramines considered in 
this study is assumed. 

In the present treatment, the intermolecular interactions between the molecules of the crystal 
as the sum of pairwise Buckingham (6-exp) (repulsion and dispersion) and Coulombic (C) 
potentials is approximated: 


V6-exp 
v ap 


(r)=A op exp(-B ap r)-C a p/r 6 , 


( 1 ) 


and 



4718 0 r ’ 


( 2 ) 


where r is the interatomic distance between atoms a and p, q a and qp are the electrostatic charges 
on the atoms, and so is the dielectric permittivity constant of vacuum. 


The parameters for the 6-exp potential in equation (1) are those previously determined for the 
case of the RDX crystal [1]. The same combination rules are used for the calculation of the 
heteroatom parameters from homoatom parameters, as previously reported [1], 

The assignment of the electrostatic charges is made by using the set of atom-centered 
monopole charges for the isolated molecule that best reproduces the quantum mechanically 
derived electrostatic potential, which is calculated over grid points surrounding the van der 
Waals surface of the molecules. This method of fitting the electrostatic potential was proposed 
by Breneman and Wiberg [13] and is incorporated in the Gaussian 94 package of programs [21] 
under the keyword CHELPG (electrostatic-potential-derived atomic charges). This method has 
the advantage of a higher density of points and a better selection procedure, which ensures a 
significant decrease in orientation effects compared to those observed with similar methods [12]. 
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The CHELPG charges were found to be invariant to either the rotation of the molecular 
coordinates or internal bond rotations. These calculations have been done at both the HF [14] 
and second-order MP2 [15-18] levels to investigate the effect of electron correlation. 

For the purpose of comparison and as an alternative to the computationally demanding MP2 
method, density functional theory (DFT) has also been used in the Kohn-Sham formulation [22], 
The DFT methods offer a less expensive but still accurate computational alternative to ab initio 
methods for including the electron correlation in post-HF treatments. In particular, the exchange 
functio nal, described by the fitted three-parameter hybrid of Becke [16] and the correlation 
functional of Lee, Yang, and Parr (B3LYP) [17] was employed. All the aformentioned theoretical 
calculations were done using a reasonable quality basis set (i.e., 6-31G** [split-valence plus d- 
type and p-type polarization functions]) [24]. 

It has been shown previously [25, 26] that the neglect of electron correlation in self-consistent 
wave functions overestimates the electrostatic interactions; however, this is mainly a scaling 
effect. Cox and Williams [25] have suggested that a scaling factor of 0.9 can be used to improve 
agreement between the calculated and experimental values of the dipole moments for a set of 
eight small molecules. The same factor has been justified in a study of the electrostatic 
interactions of a dipeptide [26], as well as in a more recent work related to the role of 
electrostatic interactions in determining the crystal structures of polar organic molecules [10]. 
Such an electrostatic model has been employed to further evaluate the effects of this scaling 
procedure. Specifically, four electrostatic models were tested for each of the 30 crystals. Two of 
them use electron correlation methods (namely MP2 and B3LYP); the third one uses unsealed HF 
charges, and the last HF charges scaled by 0.9 (denoted as 0.9HF). 

3. Computational Approach 

A general procedure for testing intermolecular potential energy functions for organic crystals 
is based on the use of molecular packing calculations [5, 6]. The basic idea is to minimize the 
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lattice energy with respect to the structural degrees of freedom of the crystal. For crystals with 
one molecule in the asymmetric unit occupying an arbitrary position, the maximum number of 
degrees of freedom is 12 and corresponds to the six unit cell constants (a, b, c, a, 3, y), the three 
rotations (0i, 0 2 , 0 3 ) and the three translations (xi, x 2 , x 3 ) of the rigid molecule. A reduced 
number of structural degrees of freedom might be involved, depending on the symmetry 
restrictions of different space groups. For crystals with more than one molecule in the 
asymmetric unit, additional degrees of freedom are introduced to describe the rotation and 
translation of the additional molecules. 

Assuming that the crystal energy is known as a function of the structural lattice parameters, 
the equilibrium crystal configuration is determined by the conditions of zero force and torque, 
together with the requirement that there is a minimum. The search for such a minimum can be 
done using a combination of steepest-descent and Newton-Raphson procedures [27,28]. 

In the present study, it is assumed that the crystals can be represented as an ensemble of rigid 
molecules. The minimization of the lattice energy with symmetry constraints has been 
performed using the molecular packing program PCK91 [29] by taking the experimentally 
observed geometries as starting configurations. This program employs an accelerated 
convergence method [1, 28] for accurate evaluation of the crystal Coulombic and dispersion 
lattice sums, with the first and second derivatives of the crystal energy evaluated analytically. In 
all calculations, a cutoff distance of 19 A has been used with the parameters T| that determines 
the relative contributions of the real- and reciprocal-space terms as defined in Sorescu et al. [1] 
having the values t|i = 0.1861 A 1 and T|6= 0.2304 A L The space-group symmetry is maintained 
throughout the energy minimization. This reduces the number of independent variables in the 
minimization procedure, resulting in a significant decrease in the computational time compared 
to unconstrained energy minimization. For example, for the 3-phase of HMX crystal with space 
group P2]/n (Z = 2), the crystallographic parameters that were varied in the minimization using 
the PCK91 program are the three dimensions of the unit cell and the angle 3, while the angles 
ct and y were frozen at 90°. Since the molecule in the asymmetric unit occupies an inversion 
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center, only the three rotations of this molecule were allowed to vary, while the three translations 
were not modified to maintain the symmetry imposed by the inversion center. It has been 
previously shown [1, 2, 30] that, despite the symmetry restrictions imposed in the PCK91 
program, the final lattice energies and crystallographic parameters are in good agreement with 
those obtained when symmetry constraints were removed (i.e., the predicted crystals maintained 
the observed space-group symmetry). 

The quality of the predicted geometrical crystallographic parameters relative to the 
experimental values has been done using a structural shift factor of the form [28,31] 


F = (A0/2) 2 + (10Ax) 2 + (100Aa/a) 2 + (100Ab/b) 2 + (100Ac/c) 2 + (Act) 2 + (A0) 2 + (Ay) 2 , (3) 

where A0 is the total root-mean-square (rms) rigid-body rotational displacement (in degrees) after 
minimization; Ax is the rms total rigid-body translational displacement (in angstroms); and a, b, c, 
and a, p, y are, respectively, the lengths of the edges and the angles of the unit cell. 

An important test of the validity of the model is the accuracy of the predicted lattice energies 
of the crystals. The lattice energies determine the relative stabilities of the different 
crystallographic phases. The calculated static lattice energy can be compared to the experimental 
sublimation enthalpy by the using the relationship [32] - AHsubi = E + Ko + 2RT, where E is the 
lattice energy and Ko is the zero-point energy. Often, a rough estimation of the lattice energy is 
obtained by neglecting the Ko term. Kitaigorodski [5] has pointed^ out that, considering the 
inaccuracy involved in the experimental determination of AHa* and with the neglect of zero-point 
energy, discrepancies up to 3-4 kcal/mol between the calculated and the observed enthalpies of 
sublimation are expected [5]. In the case of RDX crystals, it was found, using this approximation, 
that the predicted lattice energy (E = -130.09 kJ/mol) is in very close agreement with the 
experimental sublimation enthalphy (-AHsub = -130.1 kJ/mol) [1]. 
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4. Results and Discussions 


The 30 nitramine molecules considered in this study are shown in Figure 1. They were chosen 
as representative examples of important acyclic and cyclic nitramines. Different types of mono- 
and polycyclic nitramines have been included, particularly crystals that are important energetic 
materials. The structures of most of these crystals have been obtained by x-ray diffraction 
techniques. Despite the generally poorer resolution of hydrogen atom positions obtained by these 
techniques, no additional adjustment of these positions has been done to give, for example, the 
standard bond lengths [33]. The crystal structures in Figure 1 are denoted first by the common 
names of the molecules. Where available, the crystal abbreviation from the original reference is 
also included with the corresponding crystal “refcode” used in the Cambridge Structural Database 

[34] and indicated by the term in the second set of parentheses. The structures used for HNIW 

[35] and fi-HMX [36] crystals are not in the Cambridge database, so they do not have a refcode. 
In addition, different crystallographic phases of HMX and HNIW crystals have been studied; 
these are detailed in Table 1. The specific references for each of the 30 crystals in Table 1 are 
provided in brackets in the “Crystal” column and found in the reference section [35-57]. 

The results of MP calculations using the PCK91 program are presented in Table 1. The 
predicted structural lattice parameters for the great majority of these crystals deviate by less than 
2% from the experimental structures. Also, for the majority of the crystals, there are s mall 
rotations and practically no translations for the molecules in the asymmetric unit cell. The 
accuracy of the predictions can be seen in Figure 2, where the overall structural drift factors 
described in equation (3) are given. Only 10% of the total number of crystals considered here 
have a structural shift factor larger than 2.0, and practically half of the crystals have shift factors 
that are less than 1.0. 

It is important to point out that, ideally, the predicted lattice structural parameters should be 
compared with the values determined at zero temperature. However, this is not possible due to 
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Figure 1. Illustration of Molecules Whose Crystal Structures Were Studied. The Common 
Abbreviation of the Crystal Name Is Given in the First Set of Parentheses, and 
the Refcode Entry of the Cambridge Structural Database [34] Is Given in the 
Second Set of Parentheses. 
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Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures 



, values in parentheses represent the procentual difference relative to experimental values. 

For monoclinic systems, only angles are given. For triclinic systems, the cell angles and the cell angle errors are given in the order. For orthorhombic crystal 
systems, a value of 90.0 is indicated only for the experimental structure. 

J x = the net rotation of the molecule and the net translation of the center of mass of the molecule. 

F = structural shift factor; defined in equation (3). 



























































Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 



















































Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 



i values in parentheses represent the procentual difference relative to experimental values. 

For monoclinic systems, only angles are given. For triclinic systems, the cell angles and the cell angle errors are given in the order. For orthorhombic crystal 
systems, a value of 90.0 is indicated only for the experimental structure. 

5 x = the net rotation of the molecule and the net translation of the center of mass of the molecule. 

F = structural shift factor; defined in equation (3). 




























Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 
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Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 



a values in parentheses represent the procentual difference relative to experimental values. 

For monoclinic systems, only angles are given. For triclinic systems, the cell angles and the cell angle errors are given in the order. For orthorhombic crystal 
systems, a value of 90.0 is indicated only for the experimental structure. 
b x = the net rotation of the molecule and the net translation of the center of mass of the molecule. 

C F = structural shift factor; defined in equation (3). 


















































Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 





























































Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 
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For monoclinic systems, only angles are given. For triclinic systems, the cell angles and the cell angle errors are given in the order. For orthorhombic crystal 
systems, a value of 90.0 is indicated only for the experimental structure. 

b x = the net rotation of the molecule and the net translation of the center of mass of the molecule. 

C F = structural shift factor; defined in equation (3). 
d For P31 and P61, the angle = 120 is indicated. 





































Table 1. Comparison of the Crystallographic Parameters as Determined in Molecular Packing Calculations to the 
Corresponding Experimental Values for Different Molecular Structures (continued) 
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Crystal Index 

Figure 2. The Calculated Structural Shift Factor F (Equation [31) for the Crystal 
Structures. The Crystal Index Corresponds to the Number of the Crystal in 
Table 1. The Horizontal Lines at 1 and 2 Are Marked for a Clearer View of the 
Distribution of Points. 

lack of data at low temperatures. Consequently, the aforementioned comparison considers the 
deviations of the predicted geometrical parameters from the experimental values obtained at 
room temperature. It can be observed from the data given in Table 1 that the predicted lattice 
dimensions either underestimate or overestimate the experimental values. Consequently, there is 
not a general trend of the relationship between the predicted and the experimental geometric 
lattice periods, despite the small deviations between the two sets of values. 

The influence of the level of ab initio calculations on the final crystallographic parameters is 
also illustrated by the results in Table 1. The difference in the predicted geometrical parameters 
is less than 1.0% when correlated and uncorrelated ab initio methods are used. In addition, there 
is not a clear trend of the degree of accuracy with the ab initio level of calculations. In 19 of the 
systems, the accuracy increases when the HF charges are replaced with those obtained at the 
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B3LYP and MP2 levels, but, in the other 11 cases, the accuracy decreases. The maximum 
difference of the structural shift factors with the different levels of calculation is less than 0.27%. 

It can also be seen from the results in Table 1 that, when the electrostatic charges calculated 
at HF level are scaled by the 0.9 factor, the corresponding predicted geometrical parameters are 
very close to those obtained at MP2 level. Moreover, the corresponding structural shift factors 
have values intermediate between the MP2 and HF values. 

The lattice energies predicted by different models are given in Table 2. As can be seen by 
comparing the results for MP2, B3LYP, and HF methods, the use of the correlated methods 
result in a decrease in the absolute lattice energy. This effect can be understood as a consequence 
of the decrease in the absolute value of the electrostatic energy, which is attractive. The 
variations in the absolutes values of the HF lattice energies are between 8.5 and 17.5%, relative 
to the MP2 energies, with the average deviation 12.8% (see Figure 3). The use of the 0.9 scaling 
factor reduces these deviations to the range 0-7.8%, with the average deviation 4.1%. Finally, 
the B3LYP lattice energies are, as expected, much closer to the MP2 energies, with the range of 
variations 1.5-3.9% and average deviation 2.6%. These results indicate that the lattice energies 
differ significantly for sets of electrostatic charges calculated with ab initio methods that do not 
include electron correlation. These differences can be decreased by a factor of ~3 by scaling the 
HF charges. Another important result is that DFT can provide charges that give an accuracy 
(within 2.6%) for the lattice energy that is comparable to those determined at the MP2 level. 
These results are important since the computational times for B3LYP calculations are 
significantly lower than those for MP2. 

The calculated lattice energies can also be compared to experimental sublimation enthalpies 
in Table 2. For RDX (CTMTNA), HMX (0-HMX, OCHTET03) and dimethylnitramine 
(METNAM08) crystals the agreement of the MP2 energies to the experimental values is very 
good, while for tetryl (NTNANL) the difference of 15.7 kJ/mol is within the range 12-17 kJ/mol 
considered acceptable by Kitaigorodsky [5]. Despite the limited number of experimental values 
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Table 2. Comparison of the Experimental and Calculated Lattice Energies for Different 
Sets of Electrostatic Charges 


IncT 

Crystal 

AH sub 

(kJ/mol) 

Lattice Energies 
(kJ/mol) 

MP2/6-31G** 

B3LYP/6-31G** 

0.9 HF/6-31G** 

HF/6-31G** 

T 

CTMTNA 

130.1 [58] 

-130.09 

-133.68 

-136.05 

-148.18 

2 

CIWMEA10 


-113.01 

-116.40 

-118.81 

-129.07 

3 

KOFKAR 


-116.55 

-119.40 

-119.61 

-129.26 

4 

NOHTAZ 


-141.45 

-145.35 

-145.99 

-157.25 

5 

KOFKEV 


-120.00 

-122.66 

-122.20 

-131.95 

6 

p-HMX 

175.2 [58] 

-180.23 

-185.02 

-187.01 

-204.45 

7 

OCHTET 


-179.15 

-184.66 

-188.74 

-208.75 

8 

OCHTET03 

161.9 [59] 

-168.24 

-173.44 

-178.11 

-197.67 

9 

CATJIQ 


-182.50 

-187.59 

-196.88 

-213.26 

10 

JEXLUT 


-178.54 

-184.32 

-186.70 

-200.43 

11 

JEXMH 


-169.22 

-172.11 

-177.36 

-190.98 

12 

JEXMEE 


-172.98 

-178.48 

-179.86 

-196.13 

13 

KEMTIF 


-184.42 

-190.46 

-192.56 

-208.80 

14 

JEXMAA 


-176.33 

-182.29 

181.97 

-196.64 

15 

SECVOL 


-187.66 

-192.30 

-193.63 

-209.42 

16 

HNIW (e-phase) 


-186.77 

-192.82 

-196.52 

-210.88 

17 

HNIW (p-phase) 


-181.29 

-186.28 

-190.54 

-201.81 

18 

HNIW (y-phase) 


-175.31 

-180.90 

-186.15 

-198.61 

19 

DNPMTA 


-140.60 

-142.78 

-142.60 

-154.30 

20 

MTNANL 

133.8 ± 1.6 [60] 

-149.53 

-153.81 

-158.30 

-170.57 

21 

KOFKIZ 


-262.18 

-267.98 

-265.78 

-288.98 

22 

JEDSUG 


-124.22 

-128.89 

-131.24 

-145.14 

23 

JEHLAJ 


-180.70 

-187.76 

-192.29 

-209.63 

24 

METNAM08 

69.87 [61] 

-70.26 

-71.20 

-70.21 

-76.25 

25 

GEJXAU 


-171.87 

-174.31 

-175.69 

-191.41 

26 

NXENAM01 


-131.84 

-134.80 

-137.57 

-148.40 

27 

NABMUY01 


-169.28 

-173.70 

-180.86 

-196.26 

28 

DILFUZ 


-216.09 

-220.50 

-209.57 

-241.50 

29 

NOETNA02 


-166.55 

-163.68 

-172.29 

-181.04 

30 

ACENIH 


-186.20 

-188.32 

-191.61 

-205.65 


available for comparisons, it can be seen that a significant improvement in the accuracy of the 
predicted lattice energies can be obtained by using the electrostatic charges determined by 
electron correlated methods. The scaling of the HF charges also leads to improvements of the 
predicted energies, but the differences from the experimental values are larger than those obtained 
when the charges are calculated with electron correlation methods. 
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Figure 3. The Percent Differences Between the Lattice Energies and Those Based on the 
MP2 Values. The Crystal Index Corresponds to the Number of the Crystals in 
Table 1. The Three Horizontal Lines Indicate the Average Deviations for the 
Energies Calculated Using the B3LYP (<pl> = 2.6%), 0.9*HF (<p2> = 4.1%) 
and HF (<p3> = 12.8%) Sets of Charges. 
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The relative stability of different polymorphic phases of the HMX and HNIW crystals has also 
been investigated. The calculated MP2 lattice energies for the 0, a, and 8 phases of HMX are 
-180.23, -179.15 and -168.24 kJ/mol, respectively. These values support the polymorph 
stability ranking (3 > a > 8 found experimentally by McCrone [62]. Also, the calculated lattice 
energies per molecule for the e-, 0-, and y-HNIW phases of -186.77, -181.29, and 
-175.31 kJ/mol, respectively, are consistent with the stability ranking e > 0 > y reported by 
Russell et al. [4], 

5. Conclusions 

An investigation has been done on the transferability of a 6-exp Buckingham potential 
previously developed for the a-RDX crystal [1] to 30 crystals, consisting of acyclic, monocyclic, 
or polycyclic nitramines. The intermolecular potential includes Coulombic interactions between 
electrostatic charges. These charges were determined from fits to ab initio electrostatic potentials 
calculated for the individual molecules in the experimental configurations. 

The tests of this potential for the set of 30 crystals have been performed using molecular 
packing calculations. Accurate values of the crystal lattice energy have been obtained by 
employing the accelerated convergence technique for the dispersion and Coulombic lattice sums. 
Four different electrostatic models, have been considered with charges determined at HF, B3LYP, 
and MP2 levels, and with charges obtained at HF level uniformly scaled by a factor of 0.9. The 
predicted geometries indicate a good agreement with the experimental values for the great 
majority of the crystals in the study. For 90% of the crystals, the structural shift factor was less 
than 2.0, while for 50% of them, it is less than 1.0. 

There is only a small influence, generally below 1%, on the crystallographic parameters by the 
set of electrostatic charges used. However, the lattice energies are strongly dependent on the 
electrostatic model. In particular, the best overall agreement with the experimental lattice 
energies was obtained by using MP2 calculated charges. The lattice energies calculated using the 
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B3LYP charges overestimate the MP2 energies by about 2.6%, while the overestimation in the 
case of HF charges is about 12.8%. The procedure of uniformly scaling the HF charges, [10,25] 
decreases the differences to about 4.1%. It was also shown that this intermolecular potential 
correctly describes the order of stability of different phases. The predicted stability P > a > 5 is 
in accord with the experimental findings [62]. Also, the calculated stability ranking £ > P > y for 
HNIW agrees with the previously reported results [4]. 

The success of the present potential energy parameters in describing different types of 
nitramines and different phases at moderate temperatures and low pressure provides incentive to 
further investigate the transferability of this model to other classes of crystals. Incorporation of 
intramolecular motion by relaxing the rigid molecular model will also be investigated in future 
studies. 
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